Computing the optimal protocol for finite-time processes in stochastic 

thermodynamics 
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Asking for the optimal protocol of an external control parameter that minimizes the mean work 
required to drive a nano-scale system from one equilibrium state to another in finite time, Schmiedl 
and Seifert (Phys. Rev. Lett. 98, 108301 (2007)) found the Euler-Lagrange equation to be a non- 
local integro-differential equation of correlation functions. For two linear examples, we show how 
this integro-differential equation can be solved analytically. For non-linear physical systems we show 
how the optimal protocol can be found numerically and demonstrate that there may exist several 
distinct optimal protocols simultaneously, and we present optimal protocols that have one, two, and 
three jumps, respectively. 

PACS numbers: 05.40.-a, 82.70.Dd, 87.15. He, 05.70.Ln 

Keywords: fluctuation phenomena, colloids, dynamics, nonequilibrium thermodynamics 



in 



i 

C 

o 
o 



> 

cn 



o 



X 



I. INTRODUCTION 

There exist plenty of reasons for processes to be opti- 
mized. Both, economically and ecologically, it is of high 
interest to minimize the energy consumption. Mathe- 
matically, the issue addresses the question of designing 
an optimal protocol X(t) according to which some dy- 
namical system is driven from a given initial state to 
some other desired final state. Macroscopically, from the 
second law of thermodynamics, there is a lower bound on 
the required work. If started from thermal equilibrium, 
the applied work needed to reach the final state in an 
isothermal process is always larger than or equal to the 
free energy difference. The amount of the applied work 
that exceeds the free energy difference is called the dissi- 
pated work and is lost by heating the environment. Since 
tasks need to be finished within a given finite amount of 
time, the dissipated work is always positive and depends 
on the details of the protocol resulting in a technical chal- 
lenge of how to prevent powerful engines and fast micro- 
processors from their heat death. 

Along with miniaturization new aspects arise. Micro- 
scopic systems are subject to both, deterministic and 
stochastic forces, and it becomes necessary to consider 
ensemble averages. To ensure that a quantum computer 
is in a well defined state, the fluctuations have to be mini- 
mized. Finding the optimal protocol that yields the least 
fluctuations (beyond cooling down the system) is also of 
interest in soft and biomatter systems (where cooling is 
not even be possible) . In these situations the second law 
only yields constraints for the average behavior, whereas 
individual realizations may extract work from the heat 
bath thereby consuming instead of producing entropy [l[ . 
The characterization of work and heat distributions in 
fluctuating non-equilibrium situations has benefited from 
recent progress in statistical mechanics centered around 
the so-called work and fluctuation theorems, see 0, H, 0] 
and 

@, IE SB respectively. 
A trademark of these identities are exponential aver- 
ages which are dominated by the large deviation proper- 



ties of the underlying probability distributions. When 
implementing, e.g., the Jarzynski equation, e~@ AF — 
(e~@ w ) to estimate the free energy difference AF from 
the distribution P(W) of the work, where j3 = is 
the inverse temperature, high accuracy is hindered by 
the fact that P(W) markedly differs from the distribu- 
tion P(W) oc e~ f3w P(W) of dominant contributions to 
the exponential average [10]. It is then relevant to ask 
for the protocol X(t) that minimizes the mismatch be- 
tween P(W) and P(W). A convenient measure for the 
similarity between the two distributions is the so-called 
Kullback-Leibler divergence [ll[ 
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Consequently one has to search for the protocol that min- 
imizes the average work. This statement can be sharp- 
ened in the linear response regime, where the work fluc- 
tuations are proportional to the dissipated work 0, [l^] , 
(AW 2 ) — i(Wdiss)) as proven in [lj]. If one is inter- 
ested in extracting a sharp value for the mean work that 
is required, e.g. for folding a protein, it is desirable to 
drive the folding according to a protocol that minimizes 
the fluctuations of the work, (AW 2 ) = (W 2 ) - (W) 2 = 
■g(Wdiss)) resulting again in the protocol that minimizes 
the average work. (Far from equilibrium the connection 
between work fluctuations and dissipation is more com- 
plicated and no general and precise relation between the 
two is known.) 

In the present paper we investigate the uniqueness and 
general properties of optimal protocols A(i) which min- 
imize the average work necessary to accomplish a given 
isothermal transition between two equilibrium states. We 
first rederive the results obtained in [14| for linear sys- 
tems. We then consider a non-linear system and study 
the behavior well beyond the validity of linear response. 
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A. Model specifics 

We are seeking for the optimal protocol X(t) that mini- 
mizes the average work in a finite-time process of a small 
system that is subject to both, deterministic and stochas- 
tic forces. The former can be controlled experimentally 
by the external parameter A(t), whereas the latter origi- 
nate from thermal fluctuations of the environment. Tak- 
ing the environment to be a heat bath in thermal equilib- 
rium at temperature T, the stochastic forces are modeled 

via a Gaussian white noise, y-§C(£)j that is characterized 
by its vanishing ensemble average, (£(t)) = 0, and by the 
absence of any time-correlations, (C(*)C(0) = — 
Concerning the deterministic forces, we assume that 
Stokes friction is present, F — — jx, where the force of 
friction is proportional to the velocity x of the system 
along its trajectory x(t). The proportionality constant 
7 is the friction coefficient. All the other deterministic 
forces are assumed to be conservative. 

For systems on the nanometer scale-size, e.g. biologi- 
cal systems on the cellular or subcellular level and single 
molecule experiments, the stochastic forces exceed the 
inertial forces by far. Neglecting the acceleration term 
mx in the overdamped motion, where m is the mass that 
is accelerated, the microscopic dynamics is described by 
the Langevin equation 



(2) 



1 x = -V x V{x{t),X{t)) + JU{t) 



where V — V(x, A) is the potential of the conservative 
deterministic forces whose time-dependence is attributed 
to the control parameter A(t). Rescaling the time, we 
can set the friction constant to unity, 7=1. 

The time evolution of the probability distribution 
p(x, t) to observe the system at position x at time t is 
governed by the Fokker-Planck equation 



d tP (x,t) = v x ( P v x v) + ±v 2 xP . 



(3) 



Starting in thermal equilibrium, the initial canonical dis- 
tribution is 



p(x,0) = 



-0V{x,\ o ) 



(4) 



where the normalization constant Z is the partition 
function. According to the stochastic forces, any tra- 
jectory < t < tf, is possible and occurs with 
probability 

p[x(-), A(-)] = Afcxp(~ f S dt'L(x, x, A)), (5) 



where the integrand in the exponent of the probability 
functional reads 



and Af is a normalization constant. Each realization of 
the process requires its specific amount of work 



where we take possible jumps at the beginning and at 
the end of the protocol explicitly into account, e.g. 

v(rS f = V(x(t f ),X(t+)) - V{x(t f ), X(tj)). (8) 

Averaging W^[A(-), x(-)] over the initial distribution and 
the noisy history the average work 



7 ..dV Trt t'=t+ 
At' \ L T/ / 



(W[X(-)}) = / dx oP (x ,0) / dx f x 

x [ Xf ' tf Vx(-)p[x(-),H-WlX(-),x(-)}, (9) 

J(x o ,0) 

becomes a functional of the protocol [X(t)], < t < tf, 
according to which the parameter X(t) is varied from its 
initial value Ao at t — to its final value A/ at t = t/. 

The optimal protocol is found by solving the non-local 
Euler-Lagrange equation 



= 



(W[X(.)}), 



SX(t) 

where the variation reads [bj ] 

S d dV d 2 V 



(10) 
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B. Known results 

Schmiedl and Seifert studied the motion of a colloidal 
particle in an optical tweezer [l4| . For two cases, namely 
for varying the position and the strength of the trap, 
respectively, they could express the mean work as a local 
functional of one variable. This allowed them to find the 
optimal protocol that minimizes the mean work. As a 
surprising result Schmiedl and Seifert found the optimal 
protocol to jump at the beginning and at the end of the 
process, i.e. at t = and t = tf, whereas in between the 
optimal protocol varies smoothly. The initial jump can 
be interpreted as an immediate jump from equilibrium 
to a stationary state in order not to loose valuable time 
and the final jump allows a slower driving of the system 
at earlier times. It is worth noticing that the optimal 
protocol is unique in the two cases studied by Schmiedl 
and Seifert. Both, the uniqueness of the protocol and the 
possibility to express the mean work as a local functional 
of one variable, result from the linearity of the systems 
considered. 
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Below, we show how these results can be rederived by 
an explicit solution of the Euler-Lagrange equation (TIT))) , 
(fTTjl . For a generic non- linear problem such an analytic 
solution seems impossible. We therefore analyze a sim- 
ple non-linear system numerically and discuss the new 
features of the optimal protocol arising in this case. 



II. THE ANALYTIC SOLUTION 

In the following we demonstrate that an analytic solu- 
tion of the Euler-Lagrange equation (jTUJ) , (fTTjl is possi- 
ble if, as a necessary condition, the underlying Langevin 
equation ^ can be explicitly integrated. The main trick 
consists of combining the Euler-Lagrange equation with 
its derivatives in such a way that the integrals cancel out. 
Below, this is carried through explicitly for the two cases 
studied by Schmiedl and Seifert, i.e. for the stochastic 
motion of a colloidal particle in an optical tweezer. 



A. Case study I 

Dragging a colloidal particle through a viscous fluid by 
an optical tweezer with harmonic potential 



to be equal to 



V(x,X) = h(x - A) 2 



(12) 



where the focus of the optical tweezer is moved according 
to a protocol X(t) from A(0~) = A to X(t^ ) = A/ in a 
finite time tf, the expressions appearing in the variation 
of the average work (fTTj) can be computed, 



dL 

dX\ 



*C(t), 



(13) 



x(t) = e-*(x(0) + / dr(A(r) + JjC(rW), (14) 
Jo v 

(c(^(O> = v? e ~ (t '~ t)0 (*'-*)> ( 15 ) 

and the Euler-Lagrange equation becomes 



= -e - *A((T) - e" f / df A(f )e* 
Jo 

+ 2\{t)-c t J dt , X(t')e- t -Xitpe-^f-V (16) 

for all t G (0, tf). 

Differentiating the Euler-Lagrange equation twice it is 
almost identically reproduced. From the difference be- 
tween the Euler-Lagrange equation and its second deriva- 
tive follows a simple differential equation for the optimal 
protocol, 



A(t) = Vt€ (0,*/). 



(17) 



Inserting the solution of (|17p . X(t) = at + b, back into 
(fT6"j) and setting Aq = yields the integration constants 



a = b 



A, 



t f + 2-' 

resulting in the optimal protocol to be 



(18) 



[0 , t < 0, 

Kt) = \tjbi( t + 1 ) . «<*<*,, (19) 
[x f ,*>*/, 

in agreement with the result of Schmiedl and Seifert [l4[ • 



B. Case study II 



Varying the strength of the trap, 
V(x,X) = iAz 2 , 



(20) 



with A (CP) = Ao > and A (ft) = A/ > as boundary 
conditions, the expressions appearing in the variation of 
the average work (jlip can be computed again, 



dL fli j.t \ f \ 1 



x{t) =e- A W(a;(0) + x /| / dTe A ^({ T )), 
v Jo 

( x ^t))=e- 2A ^((x 2 (0)) + l y'dre 2 ^)). 
i(|x 2 (i)} = -A(t)(x 2 (t)) + i 



with 



and 



A(t) := / dt'X{t') 
Jo 



Using Stratonovich calculus, we have 

(at) X (t)x 2 {t')) 



(21) 

(22) 

(23) 
(24) 

(25) 

(26) 



/|i e -2A(0(_5_ + f d re 2A M+4 f dre 2A ^) 
Jo Jo 



(27) 



for t' > t > 0, and the Euler-Lagrange equation reads 

1 



(28) 



where A and B are abbreviations for 

A := X(t)e- 2A W - J'' dt'X(t')e- 2A ^"> 

-(X f -X(tJ))e- 2A M (29 ) 
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and 



III. NUMERICAL METHODS 



D 



1 



2 / dt"e 2A ( r 
'o 



(30) 



Note that C = AB - 1, C = AB + AB, C = AB + 
2AB + AB, . . . are complicated, because with A and B 
they contain integrals and exponentials of integrals. 

With A = (2A-2A 2 )e- 2A and B = 2e 2A , the functions 
AB, AB, AB, AB, ... do not contain any integral. 
Our goal is to express C by these simpler functions. 

From AC-AC = A(AB-AB)-2A 2 B and BC-BC = 
B(BA — BA) — 2B 2 A follows an ordinary differential 
equation for C, 

C = AB - 1 

_ (Ad - AC + 2A 2 B)(BC - BC + 2B 2 A) i 
~ (AB - Ab)(bA- BA) 

Remember that we are interested in extremizing the 
average work 

- p c^^(W) = o Vt€(0,t,). (32) 

Consequently, C and its derivatives have to vanish iden- 
tically. This yields 



= 



aA 3 b 3 



{AB - AB)(BA- BA) 



1 



(33) 



or 



= A(Ab) 3 + (AB - AB) 2 

= 16(A 2 - 12AAA + 8A 3 A + 16A 3 - 12A 2 A 2 ). (34) 

The latter is an ordinary differential equation for the 
optimal protocol. Using Lie symmetries (l5| it can be 
decomposed and integrated resulting in A = Ai and 
A = A2,3 with 

, / x 1 , , s a — c(l + ct) f , 

M*) = ^ and X 2t3 (t)= (1 | ct)2 j - (35) 

If plugged into the Euler-Lagrange equation l[28|). Xi(t) 
fails to describe the solution for positive A . X 2 ^(t) 
solve the Euler-Lagrange equation provided the integra- 
tion constants are a = Aq and 



-f - Xft f + yjl + 2X t f + X Xft 2 f 



tf(2 + X f t f ) 



(36) 



The optimal protocol is 

'Aq , t<0, 

m = < X °-f+ 2 Ct) , 0<t<t /) (37) 
Xf ,*>*/, 
in agreement with the result of Schmiedl and Seifert [lij . 



We have shown how the Euler-Lagrange equation can 
be solved analytically. Thereby, it was essential that 
the Langevin equation can be integrated. Often, this 
equation cannot be integrated analytically. In the latter 
case, it is not even possible to express the Euler-Lagrange 
equation as an integro-differential equation in A, because 
the correlation functions cannot be evaluated explicitly, 
and one has to resort to numerical methods. 

Numerically, we have implemented a Monte Carlo al- 
gorithm that minimizes the average work directly. For a 
given protocol the average work @ is approximated us- 
ing a finite ensemble of trajectories that are discretized 
in time, {x n (t m )}, n = 1...N and m = 0...M with 
<o = and tjvr = tf. The initial distribution {x n (0)} 
according to (U) and the noisy history {( n (t m )} for each 
trajectory of the ensemble are diced using the Ziggurat 
method [16j |. For each trajectory the Langevin equation 
@ is integrated according to the Heun-scheme [l7], [3 • 
Finally, the optimal protocol is found with the threshold 
acceptance algorithm [l9| . 

In the threshold algorithm we approximate the 
protocol by a polygon line that connects the 2Q + 2 
points {(0, A Q -), (0, A +), (±t f , X x ~), (±t f ,X 1+ ), 

A (Q-i)-)> (^cT*/' A (Q-i)+)> (*/) 
where the boundary values are A - = A and A^+ = Xf, 
and Q is a positive integer. The threshold algorithm 
is an iterative algorithm that starts with a random 
choice of initial values for {A +, . . . , Xf-} and computes 
the corresponding average work. In each iteration the 
values {A + , . . . , Ay- } are randomly perturbed and 
the corresponding average work is compared with the 
average work of the best protocol that was found in the 
previous iterations. If the protocol results in an average 
work that is not deteriorated by more than some given 
threshold value, it is used as the best protocol in the 
next iteration. Whenever no better protocol is found in 
some finite number of iterations, the algorithm lowers 
the threshold and continues iterating. If finally the 
algorithm does not find a better protocol in some finite 
number of iterations at zero threshold, it eliminates 
intermediate jumps in the protocol provided the protocol 
keeps optimal. 

It is worth to note that we dice the initial distribution 
and the noisy history only once and reuse these values 
in any iteration of the algorithm. This is not just to 
speed up the numerics, it is necessary for the algorithm 
to converge. 

Using up to about N — 7000 trajectories discretized 
into M = 154 time-steps and approximating the proto- 
col using a few points for the polygon, 2Q + 2 = 44 + 2, 
the algorithm is reasonably fast and can easily be run on 
a single processor. Starting from a random initial proto- 
col, it typically converges in less than 30000 iterations. 
Depending on the parameter values, it finds the optimal 
protocol within a few seconds or minutes, but for some 
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parameter values it runs for several hours. The result is 
a first approximation of the optimal protocol. 

In order to increase the numerical resolution, we re- 
run the algorithm using up to N = 50000 trajectories 
discretized into M = 700 time-steps. Starting from 
the previously found first approximation of the opti- 
mal protocol, but now approximating the protocol using 
2Q + 2 = 200 + 2 polygon points, the algorithm typically 
converges in less than 70000 iterations. The required 
CPU-time of the rerun is by a factor of ~ 40 larger, be- 
cause of the increased number of trajectories and time- 
steps involved. 

With the numerics at hand, we study a further exam- 
ple. 



A. Case study III 

Consider the stochastic motion of a small dipole in 
a viscous liquid driven by an external field where the 
direction of the external field is changed according to a 
protocol A(t). The protocol may start in any arbitrary 
initial direction Ao at time t — and ends in the final 
direction Xf = Ao + A A at time t — tf. For notational 
simplicity, we regard one degree of freedom and explore 
the stochastic motion 

x = ~V(x(t) t X{t)) + y/2DC{t) (38) 

in the potential 

V(x, A) = —H coa(x — A), (39) 

where the angle x characterizes the orientation of the 
dipole. In the numerics, we set the amplitude of the field 
to unity, H = 1, and stop the protocol at tf = 1. 

Being subject to two time-scales, the relaxation time 
and the time for driving the system, the motion of the 
dipole depends qualitatively on the diffusion constant, 
D = 4, which we take as an additional constant param- 
eter into account. 

Simulating the stochastic motion of the dipole, we find 
optimal protocols as displayed in figures [1] to [BJ 

Let us first discuss the numerical solutions with the 
diffusion constant kept small, D — 10 -5 . If AA = A/ — 
Ao is also small, the potential can be Taylor expanded 
around its minimum and is approximately equal to that 
of the moving laser trap of case study I. The solid line in 
figure [1] displays the optimal protocol for AA = 1.2. 

If the angle AA increases, the dynamics starts to ex- 
perience the non-linearities of the potential. For e.g. 
AA = 2.0 the final jump becomes much larger than the 
initial jump, see the solid line in figure O 

In rare cases, the numerical algorithm converges to 
another protocol, see the dotted lines in figures [T] and 
[51 The values of the average work tell us that these 
latter protocols are slightly suboptimal, (W) = 0.4620 
and (W) = 1.1792, while the optimal protocols require 




t 1 



FIG. 1: The optimal protocol (solid line) that minimizes the 
average work for V = — cos(a; — A) with D — 1CF J and AA = 
1.2 is unique and jumps twice. But it is interesting to see 
that there exists a suboptimal protocol (dotted line) whose 
average work, (W) = 0.4620, is close to the average work of 
the optimal protocol, (W) = 0.4619. For better visualization 
we have shifted the suboptimal protocol slightly in time. 
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FIG. 2: The optimal protocol (solid line) that minimizes the 
average work for V = — cos(a; — A) with D — 10~ J and AA = 
2.0 requires an average work of (W) = 1.1791. With (W) = 
1.1792, the average work of the suboptimal protocol (dotted 
line) almost approaches that of the optimal protocol. 



(W) = 0.4619 and (W) = 1.1791, respectively. However, 
we cannot really decide whether the dotted line in figure 
[2] is a suboptimal protocol or whether it is optimal, be- 
cause the differences in the average work are close to the 
resolution of our numerical procedure. 

If the angle AA is further increased, to e.g. AA = 3.0, 
the previously suboptimal protocol becomes optimal. 
Both protocols, indicated by the solid and the dashed 
lines in figure [3] yield the same value for the average work, 
(W) — 1.9802. The two optimal protocols differ by the 
sizes of their initial and final jumps, whereas their slopes 
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FIG. 3: For V = - cos(x - A) with D = 10" 5 and AA = 3.0 
there is a whole family of optimal protocols. Two of them 
jump twice (solid and dashed lines). All the other optimal 
protocols jump at least three times. One of them is displayed 
by the dotted line. The average work for any member of the 
family of optimal protocols is (W) = 1.9802. 



FIG. 5: For V = -cos(:e - A) with AA = n and D = 10~ 5 
the initial and the final jump vanish and the optimal protocol 
can degenerate to one single jump that may happen at any 
arbitrary time, < t < tf. Two members of the family of 
optimal protocols are displayed. The average work is (W) = 
1.99999. 




- 3.6 

- 7T 



- -0.7 
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FIG. 4: Two members of the family of optimal protocols are 
displayed for V = — cos(x — A) with D = 10 -5 and AA = 3.6. 
The two displayed protocols jump twice. All the other optimal 
protocols jump at least three times. The average work is 
(W) = 1.8130. 



A 




t 1 



FIG. 6: Optimal protocols that minimize the average work for 
V — — cos(x — A) with AA = n. For D = 0.5 there is a family 
of optimal protocols. The solid line mimics one member of 
this family. In contrast, for D = 0.7 the optimal protocol is 
unique and jumps twice as mimiced by the dashed line. 



are identical for < t < tf. 

Being the result of a numerical procedure, we can never 
claim in a strict mathematical sense that two protocols 
are both optimal. In practice, however, it is impossible 
to distinguish the case of two optimal protocols from that 
of two protocols with extremely close work values. 

The existence of two optimal protocols results from the 
symmetry in the potential, — cos(ir— A) = cos(7r— (x— A)). 
The optimal protocol, as indicated by the solid line in 
figure [31 jumps at t = to a stable state that pulls the 
dipole towards a new orientation, see figure [7] (left). This 
protocol needs most of its average work in the final jump. 



The other optimal protocol, displayed by the dashed line 
in figure[31 requires most of the average work in the initial 
jump. Immediately after this jump, the system is in a 
state where the dipole is pushed by the potential, see 
figure[7] (right). Because of the symmetry in the potential 
between its minimum and its maximum, the slopes of the 
two optimal protocols are the same for < t < tf. 

Beyond these two optimal protocols, for D small and 
AA close to 7r, there is a whole family of optimal pro- 
tocols. Any protocol of this family starts with an initial 
jump that brings the system to its pulled or pushed state, 
respectively. At any arbitrary time the protocol can jump 
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FIG. 7: The two different states that drive the orientation of 
the dipole. The left figure displays the stable state where the 
system is pulled by the minimum of the potential. The right 
figure displays the unstable state where the system is pushed 
by the maximum of the potential. 



again and the system switches between the pulled and the 
pushed state. Such a protocol is displayed by the dotted 
line in figure [3l Finally, the optimal protocol can jump 
repeatedly between the the pulled and the pushed state 
allowing for any number of jumps. 

For A A approaching tt, the size of the initial and/or 
the final jump, and the slope of the optimal protocols 
decrease towards zero, cf. figures [5] and [3] 

Increasing AA further such that it exceeds 7r, the ini- 
tial and/or the final jump, and the slope of the optimal 
protocols change sign, cf. figures [3] and [4] 

Due to the change in signs and because of the sym- 
metry and periodicity of the potential, we know that for 
A A = tt the initial and/or the final jump, and the slope of 
the optimal protocols vanish identically, and we conclude 
the following: 

For D small and AA close to tt, the optimal protocol 
jumps at least two times, except if AA = tt holds iden- 
tically. In the latter case, the initial and/or the final 
jump vanish and the whole protocol can degenerate to 
one single jump that may happen at any arbitrary time, 
< t < tt. Two members of the family of optimal pro- 
tocols are displayed in figure 03 

It is interesting to consider also the variances of the 
work for the different optimal protocols. While they give 
the same value for the average work, the variances of 
the work differ slightly among the family members. The 
optimal protocol that pulls the system all the time in 
the stable state, see figure [7] (left), yields the smallest 
variance for the work, whereas the optimal protocol that 
pushes the system all the time in the unstable state, fig- 
ure [7] (right), yields the largest variance for the work. 
Note that different values for the variances of work for 
protocols that yield the same average work imply that 
we are not in the linear response regime. 

Because of the limited numerical resolution, the slight 
differences in the variances of the work have to be treated 
with some care. Nevertheless, the fact that we are well 
beyond the linear response regime can clearly be demon- 
strated by checking whether a central identity of linear 
response theory, (AW 2 ) = 2D(W diss ) (HHGl, is vio- 
lated. For example, the optimal protocol for AA = 2.0, 
see figured results in 2D(W Aiss ) / (AW 2 ) = 1.76 ^ 1, and 
the optimal protocols for AA = it, see figure [5l result in 
2D(W diBS )/(AW 2 ) « 10 5 + 1. 

Yet, we have kept the diffusion constant D small. In- 



creasing the diffusion constant increases the noise, see 
the scrambled lines that mimic the optimal protocols in 
figure O In principle, we can get rid of the fluctuations 
in taking the average over larger ensembles, but we have 
decided to use only 50000 trajectories for the ensembles 
in order to keep the CPU-time reasonable. Moreover, in 
computing the optimal protocol for D = 0.7 we have re- 
duced the resolution to 2Q + 2 = 100+2 points, otherwise 
the CPU-time would exceed two weeks. 

Beyond the quantitative change in the amplitude of the 
fluctuations, there is also a qualitative change in the be- 
havior of the system. As discussed above, if the direction 
of the external field is changed by the angle AA = tt in 
time tf = 1 and the diffusion constant is small, there ex- 
ist optimal protocols that consist of only one single jump 
that may happen at any arbitrary time. The solid line 
in figure [5] shows that such protocols are optimal up to 
D = 0.5. 

If the diffusion becomes larger, the previously optimal 
protocols that were allowed to jump several times, get 
suboptimal. Instead, a different protocol becomes opti- 
mal. The latter is unique. The transition happens some- 
where below D = 0.7. The dashed line in figure [5] mimics 
the optimal protocol for D = 0.7, having two jumps, one 
at t = and one at t = tf = 1 . 

There are hence two regions in the AA-D-plane. These 
are shown qualitatively in figure [5J One region is located 
around AA « tt and D small. In this region, the optimal 
protocol occurs in a family and may jump several times. 
In the other region, i.e. for AA far away from tt or for D 
large, the optimal protocol is unique and jumps twice. 

If one comes close to the curve that separates the two 
regions in the AA-D-plane, it becomes hard to decide 
numerically which protocols are optimal and which are 
suboptimal, because the values of their average work ap- 
proach each other. Being blurred by the noisy history 
of the trajectories, it becomes quite impossible to deter- 
mine the exact curve that separates the two regions in 
the AA-D-plane. For this reason, we show in figure [8] 
only a crude approximation of this curve. 



IV. CONCLUSION 

The central subject of this article was to present a 
guide according to which the non-local Euler-Lagrange 
equation, a non-linear integro-differential equation of cor- 
relation functions, can be solved in order to find the opti- 
mal protocol of an external control parameter that min- 
imizes the average work for driving a small system from 
one given equilibrium to another in finite time. 

Studying the stochastic motion of a dipole where the 
direction of the external field is varied according to some 
protocol, we have found as a surprise that the optimal 
protocol may not be unique. For suitable parameter val- 
ues the protocol can jump several times. 

The reason for the non-uniqueness of the optimal pro- 
tocol results from a symmetry in the potential that allows 
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FIG. 8: There are two regions in the AA-D-plane. For AA 
near n and D small, the optimal protocol occurs in a family 
and may jump several times. For AA far away from -k or for 
D large, the optimal protocol is unique and jumps twice. 



the dominant trajectories to be near the unstable maxi- 
mum of the potential, figure [7] (right), as an alternative 
to the stable solution, where the trajectories are near the 
minimum of the potential, figure [7] (left). 

If the diffusion becomes large, the optimal protocol 
becomes unique. 

It is tempting to speculate whether biological systems 
on the cellular and subcellular level benefit from non- 
unique optimal protocols. Just imagine a protocol to be 
a triggering mechanism for a cellular process to be acti- 
vated. It can be optimal without the need of exact tim- 
ing, allowing the biological system to react spontaneously 
on the environment. 
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